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Abstract 

In this paper we present the complete derivation of the effective contour model for electrical 
discharges which appears as the asymptotic limit of the minimal streamer model for the propagation 
of electric discharges, when the electron diffusion is small. It consists of two integro-differential 
equations defined at the boundary of the plasma region: one for the motion and a second equation 
for the net charge density at the interface. We have computed explicit solutions with cylindrical 
symmetry and found the dispersion relation for small symmetry-breaking perturbations in the case 
of finite resistivity. We implement a numerical procedure to solve our model in general situations. 
As a result we compute the dispersion relation for the cylindrical case and compare it with the 
analytical predictions. Comparisons with experimental data for a 2-D positive streamers discharge 
are provided and predictions confirmed. 

PACS numbers: 51.50.+V, 52.80.-s 
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I. INTRODUCTION 



The appearance and propagation of ionization waves is the prelude of electrical breakdown 
of various media. In the case of a gas, the specific features of the breakdown waves are 
determined by the type of the gas, the value of the pressure, the geometry of the discharge cell 
and the value and variation rate of the voltage at the electrodes. The geometry determines 
the space distribution of the electric field and hence the dynamics of the ionization fronts. 
In the case where there in no initial ionization in the discharge gap, the ionization wave 
may originate from one or several overlapping electron avalanches. After attenuation of the 
electric field in the avalanche body, a conducting channel or streamer develops: a plasma 
region fully ionized with a positive side expanding towards the cathode and a negative region 
towards the anode. 

One of the approaches used to model the development of the avalanche-streamer transition 
and the streamer propagation is a nonlinear system of balance equations with a diffusion- 
drift approximation for the currents, together with Poisson equation [1]. Some progress in 
the understanding of the propagation mechanism has been achieved using that model. We 
can mention: the study of stationary plane ionization waves , self-similar solutions for 
ionization waves in cylindrical and spherical geometries , the effect of photoionization 
6] and a branching mechanism as the result of the instability of planar ionization fronts 
7-9]. In this hydrodynamic approximation, the fronts are subject to both stabilizing forces 
due to diffusion which tend to dampen out any disturbances, and destabilizing forces due 
to electric field which promote them. The solution of the model, even in the simplest cases, 
poses a challeng ing problem both numerical and analytical. Early numerical simulations can 



be found in 
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Recently, a contour dynamics model have been deduced in the limit of 



small electron diffusion 



121 ]. which resembles the Taylor- Melcher leaky dielectric model for 
electrolyte solutions [13|], but adapted to the context of electric (plasma) discharges. This 
contour dynamics model allows to study more general situations in two-dimensional and 
three-dimensional cases. 

The contour dynamics model consists of an interface separating a plasma region from a 
neutral gas region as it is shown in Fig. [U The separating surface has a net charge a and the 
thickness goes to zero as y/D being D the charge diffusion coefficient. The case displayed 
in the figure correspond a negative discharge, so the electric field is pointing towards the 
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FIG. 1. The schematic of the contour dynamics model. The case displayed corresponds to a 
negative streamer discharge discharge so a represents the negative surface charge density. The 
electric field points towards the plasma region in this case. 

plasma region and a is the negative charge density at the surface. The front will evolve 
following the equation 



where is the normal component of the electric field at the interface when approaching 
it from outside the plasma region, /i e the eletron mobility, D e is the electron diffusion 
coefficient, Eo is a characteristic ionization electric field and k the curvature of the interface. 
The parameter Iq is the microscopic ionization characteristic length. At the interface, the 
total negative surface charge density will change according to 



being now E~ the electric field at the interface coming from inside the plasma, g e is a 
parameter proportional to the resistivity of the electrons in the created plasma and j~ the 
current contribution of any electromotive force if present. 

Although the equations (0Q) and (j2J) are written for the case of a negative front plotted in 
Fig. [U and we will present the derivation of the model for this case, we could use in principle 
the same model for a positive front, but the electric field should be sign reversed, and a would 
represent the positive surface charge density. Although the moving carriers in the model are 




(1) 




da E; 



(2) 
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the electrons, one may think of a front made of holes, with a positive surface charge density, 
and characterized with the corresponding parameters for the mobility, diffusion and so on. 

In this paper, using the contour dynamics model, we will study cylindrical discharges 
when the plasma has finite conductivity. The dispersion curve for transversal instabilities 
will be obtained for these finite conductivity streamers. The results will be compared with 
the limiting cases of perfect conductivity, which is the Lozansky-Firsov model 14j with a 
correction due to electron diffusion, and with the case of a perfect insulator, i.e the limit 
of very small conductivity. Finally, we compare the results with an actual experiment for a 
positive streamer discharge. 

We start by introducing the model. Taking a minimal set of balance equations to describe 
in fully a deterministic manner the discharge (see for example we will derive the contour 
dynamics equations for the evolution of the interface between the plasma region and the gas 
region free of charge (or with a very small density of charge). The outline of this derivation 
has already been reported [12] but here we present it in full details. Then, we proceed 
by studying a cylindrical discharge in the case of finite conductivity, and the analytical 
limits of infinite resistivity and ideal conductivity. With the model at hand we will predict 
some features of the stability of the fronts. Numerical simulations are made to calculate 
the dispersion curves and test some of the analytical predictions. We briefly describe the 
numerical methods employed in the corresponding section. We end with an analysis of 
the results, the comparison with an experiment for a positive 2-D streamer discharge, and 
overview of possibilities that the model opens for more complicated geometries and fully 3D 
cases. 



II. THE DYNAMICAL CONTOUR MODEL 



In this section we obtain our model as a limit of a set of balance equations describing a 
streamer discharge. We will first recall the minimal description of a streamer discharge and 
some of the properties of the traveling planar fronts, and then make use of the asymptotic 
behaviour of those planar fronts in the limit of small diffusion to give a correction to the 
velocity of propagation of curved fronts. After finding the dynamics of the effective interface, 
a balance of the charge transport along the interface will be provided in order to complete 
the model. 
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A. The minimal model 



For simulating the dynamical streamer development of streamers out of a macroscopic 
initial ionization seed, in a non-attaching gas like argon or nitrogen, the model of a streamer 



discharge [15( can be simplified. As a first approach, the processes with the smaller prob- 
abilities or cross sections can be ignored. Attachment and recombination processes can be 
neglected on that basis in comparison with the ionization process for non-attaching gases. 
We also ignore photoionization processes in this work. With these considerations in mind, 
the resulting balance equations are 

dN 

= V • (/i e iV e E + D e WN e ) + UiN e , (3) 

w = ^- (4) 

where N e is the electron density, N p is the positive ion density, \x e is the electron mobility 
and D e the diffusion coefficient. The ionization coefficient can be modeled following the 
phenomenological approximation suggested by Townsend, which leads to 



(5) 



where Iq is the ionization length, and E is the characteristic impact ionization electric field. 



The fitting of experimental data can be done using those parameters [16] • Note also that it 
is assumed the positive ions do not move and /i e E is the drift velocity of electrons. Those 
are valid approximations at the initial stages of the streamers development, but it may not 
be right afterward. To close the model, we consider Gauss's law 

„ _ e(N„ - N e ) 

V-E = ^^ ^. (6) 

For convenience the equations are reduced to dimensionless form. Townsend approxima- 
tion provides physical scales and intrinsic parameters of the model if only impact ionization 
is present in the gas 3j. The units are given by the ionization length Iq, the character- 
istic impact ionization field E , and the electron mobility /i e . The velocity scale yields 
U = /i e E , and the time scale r = Iq/Uq. Typical values of these quantities for nitrogen 
at normal conditions are Iq ~ 2.3 /im, Eo ~ 200 kV/m, and fi e ~ 380cm 2 /Vs. We intro- 
duce the dimensionless variables r d = r// , td = t/r Q , the dimensionless field E d = E/E , 
the dimensionless electron and positive ion densities n e = N e /N and n p = N p /N with 



N = e So/(el ), and the dimensionless diffusion constant D = D e /(l U ). From now on, all 
the quantites will be dimensionless unless othewise stated. Note however that we will not 
write the subindex d. Just for reference, the dimensionless model reads 



= V • (n e E + D Vn e ) + n e a(|E|), (7) 

^ = n e a(|E|), (8) 

V • E = n p - n e , (9) 

a(|E|) = |E|exp(-l/|E|), (10) 



B. Planar fronts and boundary layer 

Using the minimal streamer model, we can compute traveling wave solutions in the planar 
case. We will assume that the plasma region is on the left and the front is moving toward the 
right. The traveling waves are solutions such as n e and n p decay exponentially at infinity. 
This means that we can take 

n e = Ae~ x{x - Vt \ 
n p = Be- X ^ x - Vt \ 
E= (E+ + Ce~ x{x - Vt) )k, 

asymptotically far ahead for the planar wave in the x direction, being E + the value of the 
electric field at the infinity. Introducing these expressions into the minimal model equations 
we get the relation 

DX 2 -{E + + v)X + a{\E + \) = 0, (11) 
which has real solutions if and only if 

v > -E + + 2^Da(\E + \). (12) 

All initial data decaying at infinity faster than Ae~ x * x , with A* = 1/ a/ Da(|E + |), will develop 
traveling waves with velocity v* = — E + + 2a/-Dq;(|E + |). Clearly, from the assumption that 
the plasma state is on the left, negative velocity solutions are unphysical. So in the case of 
a negative front, when E + is negative, the front will move at least with the drift velocity in 
the case that D — 0. For positive fronts, the motion will be possible only if the creation of 
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charge, given by the Townsend factor, and its diffusion can compensate the drift. A detailed 
discussion about the propagation mechanism can be found at |3j. 

If D <C 1 the profiles for n p and E will vary very little from the profiles with D = and 
n e will develop a boundary layer at the front. This boundary layer has a width of 0(yD) 
as shown in Fig. [2j The main results for the structure of the boundary layer which we are 
going to make use are 

n e = fix), (13) 



n 



p 



D / f(z) dz, (14) 



E = E + + 0(VD), (15) 

with x — ( x ~ v*t)/\/~D, and E the electric field in the x direction. The function fix), also 
appearing in (THl) . is the solution of the equation 

,>/ +2 Ja(|E + |)^ = /(/-!), (16) 



dx 2 V dx 

which becomes the solution of a Fisher equation under the additional assumptions that 
the Townsend factor a(|E|) ~ 1. So, as it is plotted in Fig. [21 the function / changes 
from constant values in a region of width yD, when imposing the two maching conditions 
/(— oo) = 1 and /(oo) = 0, thus separating the plasma region from the gas. The complete 
mathematical details can be found in and fl. 



C. The correction due to the curvature 



Next we will add the correction to the propagation velocity due to the curvature of the 
front. We take a level surface of n e representing the interface, and introduce local coordinates 
r (along the level surfaces of n e ) and v (orthogonal to the level surfaces of n e ). The schematic 
can be seen in Fig. El We scale the normal coordinate with the boundary layer thickness 
v = xv^D, and expand the Laplacian times D like 

da = -^ + Vd k ^ + d(a ± -^) + o(dI), 

where A = V 2 is the Laplacian operator, Aj_ is the transverse Laplacian and k is twice the 
mean curvature in 3-D or just the curvature in 2-D (details of this expansion can be found 
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n c = n 



FIG. 2. Derivation of the contour dynamics model. We take a surface of constant n e at the 
boundary which has an effective width of order y/~D. The local coordinates tangent and normal to 
the surface, r and v , together with a pillbox are also shown schematically. 



in [19J). We write ([7]) in local coordinates, and using fl9]), we find 



dn e dn e 

— 



dt 

d 2 n e 
dx 2 



dr 



dx 



+ n e a(|E|) + n e (n p - n e ) + 0(D). 



Finally we use ( IT6l) so that 

dn, 



(17) 



<)t ' Or VvD v " J d X 

Note that the curvature term correction will be relevant provided 1 k . Thus we 



have obtained a transport equation for the electron density with velocity 

v = -E + (2v/Da(|E|) - D«)n, 



(18) 



The level line n e which we have taken as representative of the interface evolution will move 
with a normal velocity 

v N = -E, + 2^/Da(\E\)-DK. (19) 

Notice that the level lines concentrate in a small region where n e presents a jump from 
its bulk value to zero, so most level lines follow ( Tr9~|) . The tangential component of the 
velocity will not change the geometry of the interface during its evolution, although tan- 
gential exchanges of charge affect the evolution through the dependence of on E v . The 
mathematical description of this effect will be the subject of next section. 
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D. Charge transport along the interface 



In order to describe the charge transport along the interface we trace a small "pillbox" 
T> around a portion of the interface having the top and bottom areas bigger than the lateral 
area, i.e. Ar ^> Av as we can see in Fig. [2J On the other hand T> will be big enough to 
contain the diffusive layer and so the portion where the total negative charge density n e — n p 
has significant values different of zero. 

We subtract ([7]) from (jSj) and integrate over the pillbox volume T>, assume that n e — > 
for x — v/VD 3> 1, |Vn e | — > for \x\ 3> 1 and get 



d_ 

dt 



f (n e - rip) dV = neE.Arl" + O(D^), (20) 
Jv 



where the contributions of the lateral transport of charge through the lateral surface is 
neglected in comparison with the exchange of charge in the normal direction. Note that in 
the Taylor- Melcher model this assumption is also made. As explained in [3] the left hand 
side of equation (1201) can be written as the time partial derivative of the product of the 
negative surface charge density a times the normal area Ar, and the change of a surface 
element can be related to the curvature times the normal velocity, so that 

da 

— + kv n o = - rieEvl^^ (21) 

If a charge source I(t) is present in the plasma, for instance at xq, this source will create 
a current density inside the plasma and we will have at the interior of Q 

V-j = /(i)5(x-x ). (22) 

By adding this contribution to (T2"T|) we can finally write 

da E~ ._ 

— + kv n u = j v , (23) 

at q 

where j~ is the current density coming from the ionized region Q to its boundary dQ in the 
normal direction u, E~ is the normal component of the electric field when approaching the 
interface from inside, and g -1 = lim x= _ 00 n e is the effective movility of the electrons inside 
the plasma. Note that the quasineutrality of the plasma, further away of the interface is not 
changed by the current, but there is a jump in the normal component of the electric field 
across the interface given by 

E+ - E; = -a, (24) 



with the normal component of the electric field when approaching the interface from 
outside the plasma region. 



E. The effective contour model 

The Eqs. (1T9]) and (1251) together constitute the dynamical model able to describe the 
evolution of an interface separating a plasma region from a neutral region. Notice that 
in the case g^ 1 3> 1, we arrive to Lozansky-Firsov model [14] with a correction due to 
electron diffusion, meanwhile in the limit D = we arrive at the classical Hele-Shaw model. 
Such a model is known to possess solutions that develop singularities in the form of cusps 
in finite time 17| but, when regularized by surface tension corrections, the interface may 
develop various patterns including some of fractal-type (see [18| for a recent development 
and references therein). 

Eq. ( 123]) will provide the surface charge density a as a function of time. From it, we 
can compute the electric field and move the interface with (119]) . Two limits can be easily 
identified in the case that there is no charge injection inside the plasma, i.e j~ ~ 0: a) the 
limit of large conductivity 

fT 1 > l, E; = o, 

so that the interface is equipotential and b) the limit of small conductivity 

d 

o -1 < 1, — (a At) = cxAr = Cte, 
at 

where the charge contained by a surface element is constant and the density only changes 
through deformation (with change of area) of the interface. In the next sections we will 
study the intermediate case of finite resistivity. 



III. THE CASE OF FINITE RESISTIVITY IN 2-D GEOMETRIES 

As an application we will solve the 2-D case for different conductivities. In order to grasp 
some features of the model first we will consider how fronts with radial symmetry evolve. 
Then we will study the stability of those fronts under small perturbations and finally solve 
the model numerically in order to test some of the analytical predictions. 
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A. Solutions with radial symmetry 

The electric potential created by a surface charge distribution with radial symmetry at 
the distance r is found by solving the equation 

AV = a8(r). (25) 

The fundamental solution turns out to be in polar coordinates 

Clog Ixl, Ixl > r 
V-(x) { (26) 

CTogr, |x| < r 

where C will be determined by the condition of the electric field jump (|24p at the surface. 
From the potential solution we can compute the electric field which has a discontinuity at 
the surface 

E~ = 0, Et = --. (27) 

r 

For the current density, the solution of fT22l) gives 

j = (28) 

and finally using ( |23l) and the fact that Vn = dr/dt and k — 1/r, we get 

da 1 dr I(t) 



This equation can be easily solved. We can write it as 

d(ra) I(t) 
dt ~ ~^T' 

to get 



(30) 



<t = -^&, with Q(t)= I I(t)dt, (31) 



2vrr Jo 
where we have assumed that <r(0) = 0. Now we can see from the condition ( 12^1) that 
C = -Q(t)/2ir, so 

E t = (32) 
Then, defining e = D, the interface evolves according to ([!]) as 



dr_ _ (Q{t) 
dt ~ 



(^- + e^)l + 2^ea(\Q(t)/2nr\). (33) 
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We shall analyze next two limiting cases. First the case where 



r« ; and e « 1. (34) 

47reV a (l<5(*)/2^l) 



Then expression (j33j) results 

_ Q(t) 

dt ~ 2vrr 

so 



(35) 



r(t) « yr(0) 2 - jT Q(f)/irdtf. (36) 
For the particular case — Q is constant 

r(t) « v/ r (0) 2 - tg/vr (37) 



The second case is the opposite one. If 



r > Jgglj _ and e<1 ( 38 ) 

47reV«(IQW/2^|) 

we have now 

^« 2e V«(I0(*)/2^D- (39) 

For the particular case = Q, by standard asymptotic calculations, when t > 1 we 
deduce 

r (t) ~M l g t (40) 

7T 

B. Stability analysis 

We will study now the stability of the fronts under small perturbations. We change by a 
small amount the position of the front as well as the charge density. The perturbed position 
and charge surface density of the interface on the interface will be parametrized using the 
polar angle as 

r(9,t)=r(t) + 6S(9,t), (41) 

^ = -d|!t) + ™ (42) 

where r(t) is the solution of the equations for the radial symmetrical front, Q(t) = J * I(t) dt 
and 5 a small parameter. 
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The electric potential will change by <5VJ,(x) after adding a geometrical perturbation of 
the interface and some extra charge on it. This term satisfies the equation AV P = 0(5). 
Changing coordinates to 

- r(t) 

x — > X = X 



the perturbed surface becomes a disk of radius r(t) again, and solving for it yieds 



oo 

'T 



V P (r, 6) = J2^n cos(n0) J , f > r (43) 
1 

V p (f, 9) = Y^pn cos(n6) (-J , f < r (44) 

where it is imposed that V p remains finite at the origin and at very large distances becomes 
zero. 

Now taking the condition of continuity for the potential, we have at the interface x s (in 
the original coordinate system) 

and writing the surface perturbation as 

oo 

S = $>„(*) cos(n0), (45) 

n=l 



the coefficients of the series in ( 143|) and (J44J) can be related by 



Q(t) , . 

Wn = ¥n + ~Z S n . (46) 



Making use of the expressions (1431) - (1461) . one can calculate the electric field to 6 order. We 
will need the normal components of the electric field at both sides of the surface, together 
with the jump condition ff24l) to find the charge perturbation of fj42|) . The normal components 
of the electric field at the interface are 

Qjt) Q(t) \n 

oo 

E^ = -6j2<Pn-cos(ne), (47) 



r 
l 



thus 



n=l v ' 



g(t) - ^ -cos(n0). (48) 
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The dynamics of the front will be changed by the perturbation introduced. The curvature 
correction turns out to be 

r 2 + 2rS5 -rS 0e 5 + O(5 2 ) 1 S + S e 



K 



-5 + 0{5 2 



(49) 



(r 2 + 2rS e 6 + 0(5 2 ))^ r r 
(the subindex 9 means the partial derivative with respect this variable) and the normal 
component of the velocity 

(50) 



dt ' dt 

so the contour model equation ffl9l) . to first order gives 



dt 



dt 



2nr 



2vrr 2 



Q(t) \n 

- — s n - cos(nfc 

ZTxr I r 



2e 2 a/«o + Sotr 



T £ S + S e 
e I d 



(51) 

where we have written the Townsend function (|TU|) up to first order as a = a + <5ai + 0(<5 2 ). 
Now, we have 

\Eq + 5Ei|e~i B o+Wi « |E |e~^T + <5sign(E )Ei(l + — )e~^T = a Q + 5a u 

\ En 



where, using ( 14T1) . 



27rr 

oo 



so that 



E^^ln^ + ^-l)-— s n l- 

n=l ^ ' 



cos(n0), 



V 7 " = a/«0 + ^ 



Ql 



2JoT n 



l + 5sign(Q(t))^fl+ 1 



2 -En 



-En 



Taking into account ( 133|) for the zero order term, we get from (!5T|) 



gg _ ? g(t) 

^ 2vrr 2 



E, . Q(t) \ n ( S + S e9 \ I o i 
( <^ n + -^S n ) - COs(nS) + £ ( 5 ) + £2 



27rr 



and finally making use of the expansion ( |4"5"1) for the perturbation S 1 yields 





i 27rr A /a^sign(Q(t)) / 2vrr \ 

■ ) + : — ^ — + 



IQWI 



n 



(52) 
(53) 



(54) 



(55) 
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In order to find the correction to the charge density we take Eq. (!23|) and multiply it by 
r(9,t). Then we use the curvature expansion (H91) written as 

1 Sqq 

K ~ KM) ~ T 2 " ' 

(being r — r(t) is the zero order term in the position), and the fact that 

dr(9,t) 



v N 



dt 



Hence 



d(r(9,t)a(9,t)) 



dt 



r(9,t)^v N a(9,t)5 



r(9,t)„_ I(t) 



-E 



2tt 



so that, at 0(5), 



0(r£) S ee Q{t)dr _ 1 



dt 



r 2nr dt g 



cosint 



Making use of the (p|) . (j45l) and (pD, we get 



d / Q(t) 
-— 2ny? n + n- — s, 
dt \ 2 r nr 



Q(t) dr 2 n 
— 2Tf n s n + -pn, 



or after simplifying 



dt 2ixr dt 
Finally, using ( )33|) and ( 1551) 



rf^n , Q(t) ds n Q[t) dr 
2— - — h — — = —- — - — (n — l)Sr 



2nr 2 dt 



I(t) 1 

S n tp r , 

lixr Q 



dt 2vrr 



n 1 27rr v / ao sign(<5(t)) / 2-nr \n 

+ «> ^ (1 + m ) -fn- 



r 



Qit) 

2nr 2 



n - l)s n -(n 2 - l)s n - e 



m)\ 

i (n - l)*/a^ 



2?rr \ 

+ JoW\)' 



Q(t)f Q(t) e i V J(t) 1 

2nr A \ zur r I 2ixr 



Q 



and after rearranging the terms 



dif n _ 1 

dt ~ 2 



— — n-£2^ 

2rcr z r 



1 + 



f g(j) 
1 2?rr 2 



2nr 



+ 



2r 



2vrr ^ 




\Q{t)\j 











l fn + 
1 



[n-V 



(56) 



(57) 



Thus the time evolution of each particular mode has been obtained and it is governed by 
fl55D and flSU). 
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C. Special limits 



First we study the limit of ideal conductivity. It corresponds to g — > 0, and hence, from 
(l5Tj) . we can conclude that tp n — > 0. Physically this means that in the limit of very high 
conductivity, the electric field inside goes to zero (E~ — > 0), as we approach to the behavior 
of a perfect conductor. If we consider that Q(t) = Qo is constant or its variation in time 
is small compared with the evolution of the modes (which also implies I(t) — > 0), and the 
same for the radius of the front r(t) = r , we can try a solution s n = exp(u n t), (p n = 0, to 
(1551) . and get a discrete dispersion relation of the form 



u n = 

2nr, 



Next we consider the limit of finite resistivity, but such that the total charge is constant 
at the surface, or varies very slowly. Writing ( )57|) as 

dip n d ( Q(t) \ Q(t) dr 



dt dt\^r Sn ) A^dt nSn 2*/ n ' 
we have now 

di£n Qo ds n ^ J_ 

dt ~ 47rr dt 2^"' 1 ' 

For a small enough conductivity, g — > oo so no extra charge reaches the surface, we find 
fn = ~£^s n , and with s n = exp(w n t), Q55J yields 

i 



^ = -^f--l)-4(- 2 -l)-^fl + ^ (--1). (60) 
2nr 2 \2 J r% } r \ \Q \ J \ 2 1 

In a curved geometry we can see that the modes are discrete. However, if we compare 

f )58|) and fl60|) . for small n and vanishingly small ao there is a 1/2 factor discrepancy in 

the dispersion curve between the two limits. The origin of this prefactor was discussed for 



planar fronts in 20(, and the dispersion relation for planar fronts was obtained in the case 
of constant charge in We get in this 2-D curved case the same factor 1/2 that we got 
for the planar case. On the other hand, imposing constant potential at the surface gives a 
factor of 1. The intermediate situations can be studied by solving the system fl55|) and floTI) . 

Another important consequence is that in both cases the maximum growth correspond 
to a perturbation with 

n^\Q \/D, (61) 
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provided that the term can be neglected and Qq is negative, implying that the number 
of fingers increases with the net charge and decreases with electron diffusion. 



D. Numerical simulations 

In order to test the analytical predictions, we have calculated numerically the dispersion 
relation curves for the cases studied previously, when g — > 0, so we have a perfect conducting 
plasma, and when q remains finite. We will outline the numerical algorithms and present 
here the results. 

We start for the case of finite resistivity. The 2-D solution for the potential problem can 
be written as 

$(x)= / -Llog|x-xV(xO^', 
Jan ^ 

Note that the integration domain dfl is the curve manifold. The electric field results 

E = -[ ±-^La(x')ds'. 
Jan ^ |x - x'| 

In order to obtain the component in the normal direction E u , we will multiply it by the 
normal pointing outside the plasma region, i.e. 



n 



x 



where the subindex denotes the derivative respect to the curve parameter /3. So we can 
write 

1 (x — x',y — y') (yfa-xp) 



Jan 2tt (x - x') 2 + (y - y') 2 J~2 + y 2 V p p 



Now when approximating the integral as a discrete sum on the interface, i.e. E+ limit, some 
care must be taken. We need the limit E+ on the interface. When x coincides with x' 
there is an extra contribution of half the pole, which is cr(x)/2. The E~ can be obtained 
from the boundary condition (j24p . and the curvature must be expressed in the appropriate 
coordinates system. 

The case of constant potential, which corresponds to g = 0, is treated numerically as 
follows. We have to fulfill the condition 



/ — log |x - x'| a(x!)ds = V Q , 
Jan 2vr 
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being Vq a constant for any x belonging to dVt. Discretizing the domain in small segments 
Ai between points x, and Xj+i we can approximate the integral as 



f 1 i 

/ 2~ lo § |Xj - X'| + — log |x"; - Xy| tr(Xj) |x i+1 - Xj| = Vb- 



3 

where Xj is the mean point of the A4 segment. The self contribution of the segment to the 
integral is taken as 

f 1 1 ,- „ , ft 1 1 , , , 1 , /. hi 
/ — log x< — x las = / — log mm = — fk log 1 

J Ai 2tt 5 1 1 2vr & 1 1 2vr V 2 

being /ij the length of Ai. So we end with the equation 

M i3 a 3 = V 1, 
where 1 is the identity matrix, Oj = cr(Xj), and 



Mi, 



±hi (log f - 1) , for i = j 



271 

Due to the linearity of the problem, we can solve A i3 a 3 = 1 and rescale subsequently the 
solution in order to fulfill a jhj = Q- 

In the numerical simulations presented here, we will follow the evolution of a total initial 
dimensionless charge Q = — 10 distributed uniformly along the curve given by 

x{9) = [1 + 0.05 cos(n#)] cos(0), 

y (6) = [1 + 0.05 cos(n6)] sin(fl). (62) 

where n gives the mode of the perturbation and 6 is the curve parameter. We assume that 
there is not input current, so j~ = in ( 123]) . and then compute the exponential growth of 
each mode for a small period of time in order to get the dispersion curve. In Fig. [31 for 
different values of the inverse of the resistivity coefficient 1/g (or effective conductivity), we 
plot the corresponding dispersion curves. 

Note that the slope increases with the increase of the conductivity of the plasma, the 
maxima moves to higher modes, and for larger n's the dispersion curves become negative as 
predicted by ( 158]) and (160]) . The slope around the origin n = 1 is larger for the case of ideal 
conductivity, i.e. when the interface is equipotential. 
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FIG. 3. Dispersion relation for the discrete modes of a perturbation with initial value Q = — 10 for 
different inverse resistivities 1/g. The ▲ are for 0, + for 5, * for 10, ■ for 15, T for 25. The case 
of zero resistivity corresponds to 

IV. COMPARISONS WITH 2-D POSITIVE DISCHARGE EXPERIMENT 



In this section we will make some estimations in order to test the validity of the assump- 
tions made in our contour dynamical model. We will use the experimental data presented at 



reference [2lJ . The experiment reported there consists in the measuring of the potential and 
electric field distribution of a surface streamer discharge on a dielectric material. For that, a 
technique based on Pockels crystals have been applied in order to obtain some temporal and 
spatial resolution of the discharge (see the reference for details). However, a note of warning 
must be done: a surface discharge is not a truly 2-D discharge, due to the fact that there 
is a vertical contribution of the electric field, and the discharge has two different interfaces, 
the air and the substrate, so the boundary conditions are not the same that the presented 
so far in this paper. Nevertheless, and keeping that in mind, we may try a quantitative 
estimation a la Fermi from our results and compare it with the actual experiment. 

Here it is a brief account of the experiment. A discharge is created on a dielectric surface 
using a positive tip and branching is observed. Then the potential is measure using Pockels 
crystals, laser pulses and a ccd camera. The temporal resolution is 3.2 ns and the electric 
field close to the tip reaches values of 3 kV/mm, leaving behind a potential gradient of 0.5 
kV/mm. At position r = 8 mm the front moves with an estimated velocity from the pictures 
of 0.18 mm/ns (from the charge density data, the front has a radius of 4 mm at 3 ns, 8 mm 
at 15 ns and 9 mm at 28 ns). The pictures show a sharp interface for the charge distribution, 
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so our model should be able to give some quantitative predictions. Unfortunately, there is 

only one discharge reported, so the estimations we are going to make are very rough. 

The experimental data gives a characteristic front speed Uq pa 0.1 mm/ns, and E ~ 200 

kV/m. In order to get an estimation of the diffusion coefficient D we can make use of the 

expression (|T9"|) . We take 

, 3kV/mm 0.185 mm/ns 
pa pa 1.5, and fjv ~ 77 ~ 1-9, 

so that .D pa 0.05 is the number that we get. Note that from the expression Uq = /i e E , we 
could find the experimantal value for the mobility /i e for this discharge. 

Now we can make a prediction. The maximum of the dispersion relation will tell us the 
number of fingers one may find in such experiments. We have calculated the dispersion 
relation for two limit cases. The limit of ideal conductivity (I58p and the limit of infinite 
resistivity (!60|) . Those limits would give a lower and upper estimation values for the actual 
dispersion relation. We expect that the experiment will lie in between and be closer to the 
predictions given by the limit of infinite resistivity, as the discharge is on a dielectric plate. 
But before using those dispersion relations we need a further estimation for the surface 
charge density. We can get the surface density from the jump of the electric field across the 
interface. So the dimensionless expression reads 

(3 - 0.5) kV/mm 

do ~ — , at r = 8mm. 

Eo 

In the dispersion relation expressions (1581) and fl60l we have to make the substitution 
<5o/2vrro = <Jq and find the maximum for n. For the ideal conductivity case (|58|) yields 
a maximum at n ~ 76, and for the case of infinite resistivity ( 160]) . turns out n pa 14. Count- 
ing the numbers of real fingers in the experimental pictures at 15 ns, the number is around 
20 (one has to extrapolate the number as the pictures do not show the whole discharge). 
This number is much closer to the lower limit as we pointed before, pointing in the direction 
that the electrons on the dielectric surface, when moving through the plasma, feel a much 
higher resistivity than in a conductor. 

Although we do not expect to capture the whole physics of the discharge with the contour 
model, some essentials ingredients for the early development of the front seem to be well 
accounted by it. The theoretical prediction made in this section is a rather good one, despite 
all the approximations made and gives some insight about the parameters involved, such 
the mobility of the carriers, diffusion coefficient, number of fingers, and so no. 
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V. CONCLUSIONS 



We have presented the complete derivation of the contour dynamics moder electric dis- 
charges introduced in [2]. The model appears as the leading asymptotic description for the 
minimal streamer model when the electron diffusion coefficient is very small. It consists of 
two integro-differential equations defined at the boundary of the plasma region: one for the 
motion of the points of the boundary where the velocity in the normal direction is given in 
terms of the electric field created by the net charge there, and a second equation for the 
evolution of the charge density at the boundary. This second equation is very similar to 
the Taylor- Melcher model in electrohydrodynamics 13). In the model the electric field is 
determined by solving Poisson equation with a given surface charge density, leading to a 
singular integral of the density. 

Once our model has been deduced, we have computed explicit solutions with cylindrical 
symmetry and investigated their stabilities. The resulting dispersion relation is such as 
the perturbation with the small mode number can grow exponentially fast. In fact, both 
the number of modes become unstable and the mode that becomes most unstable (the one 
corresponding to the dispersion relation) depends critically on the electric resistivity of the 
media. We have computed analytically the dispersion relation and found that the number 
of unstable modes grows with the inverse of the resistivity (the conductivity) and the most 
unstable mode also increases with it. In the limit of vanishing resistivity one can consider 
the medium as a perfect conductor and therefore impose that the potential is constant at 
the boundary. The dispersion relation for the model with finite resistivity converges to this 
limit when resistivity tends to zero. 

We have implemented a numerical procedure to solve our model in general situations. 
In order to develop the numerical method, we needed to evaluate certain singular integrals 
that appear when computing the electrical field. As one result of the numerical method, the 
dispersion relations have been computed and c omp ared them with the analytical results. 



As a difference to our previous communication 12], we have paid special attention to the 
effects of Townsend expression for impact ionization ([5]) on the dispersion relation and the 
cases of intermediate resistivities. 

Finally, we have taken some experimental data from a positive surface streamer discharge 
and compare them with our model predictions. The number of fingers calculated from our 
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model is of the same order of the observed one in the actual experiment. We have been able 
also to estimate the diffusion coefficient from the data. We have shown that the behaviour 
of the carriers inside the plasma is closer to the limit of high resistivity, so the importance of 
taking into account the plasma resistivity is made clear. Thus, it is proved that our contour 
model is able to capture essential parts of the physics involved in the earlier development 
the streamer discharge, with an extra bonus: we can study more complex geometries and 
general situations both analytically and numerically. We are now in the process to complete 
the fully 3-D case and extend these results. 

The authors thank support from the Spanish Ministerio de Education y Ciencia under 
projects AYA2009-14027-C05-04 and MTM2008-0325. 
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